<html>
  <head>
    <meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
    <title>derivative</title>
  </head>
  <body bgcolor="#FFFFFF">
    <center></center>
    <div align="right">Last update : 07/02/2006</div>
    <p>
      <b>derivative</b> -  approximate derivatives of a function</p>
    <h3>
      <font color="blue">Calling Sequence</font>
    </h3>
    <dl>
      <dd>
        <tt>derivative(F,x) </tt>
      </dd>
      <dd>
        <tt>[J [,H]] = derivative(F,x [,h ,order ,H_form ,Q]) </tt>
      </dd>
    </dl>
    <h3>
      <font color="blue">Parameters</font>
    </h3>
    <ul>
      <li>
        <tt>
          <b>F</b>
        </tt>: a Scilab function F: <tt>
          <b>R^n --&gt; R^m</b>
        </tt> or a <tt>
          <b>list(F,p1,...,pk)</b>
        </tt>, where F is a scilab function 
      in the form <tt>
          <b>y=F(x,p1,...,pk)</b>
        </tt>, p1, ..., pk being any scilab objects (matrices, lists,...).</li>
      <li>
        <tt>
          <b>x</b>
        </tt>: real column vector of dimension n.</li>
      <li>
        <tt>
          <b>h</b>
        </tt>: (optional) real, the stepsize used in the finite difference approximations.</li>
      <li>
        <tt>
          <b>order</b>
        </tt>: (optional) integer, the order of the finite difference formula used to approximate 
      the derivatives (order = 1,2 or 4, default is order=2 ).</li>
      <li>
        <tt>
          <b>H_form</b>
        </tt>: (optional) string, the form in which the Hessean will be returned. Possible forms are:<ul>
          <li>
            <tt>
              <b>H_form='default'  </b>
            </tt>: H is a m x (<tt>
              <b>n^2</b>
            </tt>) matrix ; in this form, the k-th row of H  corresponds to the Hessean of 
      the k-th component of F, given  as the following row vector :<pre>

  [ d(grad(F_k))/dx_1 | ..... | d(grad(F_k))/dx_n ] 
   
                </pre>
            <p>  ((grad(F_k) being a row vector).</p>
          </li>
          <li>
            <tt>
              <b>H_form='blockmat' :  </b>
            </tt>H is a (mxn) x n block matrix : the classic Hessean matrices (of each component of F) are 
    stacked by row (H = [H1 ; H2 ; ... ; Hm] in scilab syntax).</li>
          <li>
            <tt>
              <b>H_form='hypermat' :  </b>
            </tt>H is a n x n matrix for m=1, and a n x n x m hypermatrix otherwise.  H(:,:,k) is the 
    classic Hessean matrix of the k-th component of F.</li>
        </ul>
      </li>
      <li>
        <tt>
          <b>Q</b>
        </tt>: (optional) real matrix, orthogonal (default is eye(n,n)).</li>
    </ul>
    <h3>
      <font color="blue">Description</font>
    </h3>
    <p>
    Numerical approximation of the first and second derivatives of a function
    F: <tt>
        <b> R^n --&gt; R^m</b>
      </tt> at the point x. The Jacobian is computed  by approximating
    the directional derivatives of the components of F in the direction of 
    the columns of Q. 
    (For m=1, v=Q(:,k) : grad(F(x))*v = Dv(F(x)).)
    The second derivatives are computed by composition of first order derivatives.
    If H is given in its default form the Taylor series of F(x) up to terms of second 
    order is given by :</p>
    <pre>

  F(x+dx) = F(x) + J(x)*dx + 1/2*H(x)*(dx .*. dx) + ...
   
    </pre>
    <p>
    (([J,H]=derivative(F,x,H_form='default'), J=J(x), H=H(x).)</p>
    <h3>
      <font color="blue">Remarks</font>
    </h3>
    <dl>
      <p>
    Numerical approximation of derivatives is generally an unstable process.  
The step size h must be small to get a low error but if it is too small floating  
point errors will dominate by cancellation. As a rule of thumb don't change 
the  default step size.  To work around numerical difficulties one may also change 
the order and/or choose different orthogonal matrices Q (the default is eye(n,n)), 
especially if the  approximate derivatives are used in optimization routines. All 
the optional arguments may also be passed as named arguments, so that one can use 
calls in the form :</p>
      <pre>

derivative(F, x, H_form = "hypermat")
derivative(F, x, order = 4) etc.

   
    </pre>
    </dl>
    <h3>
      <font color="blue">Examples</font>
    </h3>
    <pre>

  function y=F(x)
   y=[sin(x(1)*x(2))+exp(x(2)*x(3)+x(1)) ; sum(x.^3)];
 endfunction
 function y=G(x,p) 
   y=[sin(x(1)*x(2)*p)+exp(x(2)*x(3)+x(1)) ; sum(x.^3)];
 endfunction

 x=[1;2;3];[J,H]=derivative(F,x,H_form='blockmat')

 n=3;
 // form an orthogonal matrix :   
 nu=0; while nu&lt;n, [Q,nu]=colcomp(rand(n,n)); end  
 for i=[1,2,4]
    [J,H]=derivative(F,x,order=i,H_form='blockmat',Q=Q);
    mprintf("order= %d \n",i);
    H,
 end

 p=1;h=1e-3;
 [J,H]=derivative(list(G,p),x,h,2,H_form='hypermat');H
 [J,H]=derivative(list(G,p),x,h,4,Q=Q);H

 // Taylor series example:
 dx=1e-3*[1;1;-1];
 [J,H]=derivative(F,x);
 F(x+dx)
 F(x+dx)-F(x)
 F(x+dx)-F(x)-J*dx
 F(x+dx)-F(x)-J*dx-1/2*H*(dx .*. dx)

 // A trivial example
 function y=f(x,A,p,w), y=x'*A*x+p'*x+w; endfunction
 // with Jacobian and Hessean given by J(x)=x'*(A+A')+p', and H(x)=A+A'.
 A = rand(3,3); p = rand(3,1); w = 1;
 x = rand(3,1);
 [J,H]=derivative(list(f,A,p,w),x,h=1,H_form='blockmat')
 // Since f(x) is quadratic in x, approximate derivatives of order=2 or 4 by finite
 // differences should be exact for all h~=0. The apparent errors are caused by
 // cancellation in the floating point operations, so a "big" h is choosen.
 // Comparison with the exact matrices:
 Je = x'*(A+A')+p'
 He = A+A'
 clean(Je - J)
 clean(He - H)
   
  </pre>
    <h3>
      <font color="blue">See Also</font>
    </h3>
    <p>
      <a href="numdiff.htm">
        <tt>
          <b>numdiff</b>
        </tt>
      </a>,&nbsp;&nbsp;<a href="../polynomials/derivat.htm">
        <tt>
          <b>derivat</b>
        </tt>
      </a>,&nbsp;&nbsp;</p>
    <h3>
      <font color="blue">Author</font>
    </h3>
    <p> Rainer von Seggern, Bruno Pincon</p>
  </body>
</html>
